Fast, numerically stable computation of oscillatory integrals with stationary points

We present a numerically stable way to compute oscillatory integrals of the form ∫ 1 −1 f(x)e iωg(x) dx. For each additional frequency, only a small, well-conditioned linear system with a Hessenberg matrix must be solved, and the amount of work needed decreases as the frequency increases. Moreover, we can modify the method for computing oscillatory integrals with stationary points. This is the first stable algorithm for oscillatory integrals with stationary points which does not lose accuracy as the frequency increases and does not require deformation into the complex plane.